1. Field of the Invention
The invention relates generally to 3D image reconstruction in a cone beam x-ray imaging system, and more specifically to correction of image inaccuracies caused by mask boundaries when using a simplified 3D backprojection image reconstruction technique.
2. Description of the Prior Art
A filtered backprojection (FBP) cone beam image reconstruction technique is described by Kudo, H. and Saito, T., in their article entitled "Derivation and Implementation of a Cone-Beam Reconstruction Algorithm for Nonplanar Orbits", IEEE Trans.Med. Imag., MI-13 (1994) 196-211, incorporated herein by reference.
Briefly, the FBP technique consists of the following steps at each cone beam view (i.e., at each position of the radiation source as it scans about the object, and at which an imaging detector acquires a corresponding set of projection data):
1. Compute a ID projection (i.e., line integral) of the measured cone beam image acquired on a detector plane 1 at each of a plurality of angles .theta.. This step is illustrated by FIG. 1A for a given angle .theta..sub.1 of a plurality of angles .theta., where the projection 2 at coordinates (r,.theta.) comprises the integrated values of the cone beam image 1 on detector plane 1 along plurality of parallel lines L(r,.theta.) that are normal to angle .theta., each line L being at an incremental distance r from an origin O. Generally, if the detector plane 1 comprises an N by N array of pixels, then the number of angles .theta. is typically given by .pi.N/2. PA1 2. Filter each ID projection in accordance with a d/dr filter, resulting in a new set of values at each of the r,.theta. coordinates, such as shown by filtered projection 6 for the angle .theta..sub.1 in FIG. 1A. PA1 3. Normalize the filtered projections with a normalization function M(r,.theta.). Normalization is needed to take into account the number of times the plane of integration Q(r,.theta.) which intersects the source position and the line L(r,.theta.), intersects the scan path, since the data developed at each scan path intersection creates a contribution to the image reconstruction on the plane Q(r,.theta.). PA1 4. Backproject the filtered projection 6 from each angle .theta. into a 2D object space 7 which coincides with the detector plane 1. This step is illustrated by FIG. 1B, wherein lines 8 spread the value from each r,.theta. coordinate into 2D space 7 in a direction normal to each .theta.. PA1 5. Perform a 1D d/dt filtering of the backprojection image formed in 2D space 7 by step 4. The 1D filtering is performed in the direction of the scan path, i.e., along line 10, where t points in the direction of the scan path. PA1 6. Perform a weighted 3D backprojection of the resulting data in 2D space 7 (i.e., from each pixel in the detector) onto a plurality of sample points P in a 3D object volume 12. The density assigned to each point P is weighted by the inverse of the square of the distance between the point and the spatial coordinates of the x-ray source (see Equation (59) of the forenoted Kudo et al article). PA1 1) applying a mask to each set of cone beam projection data acquired at each source position, thereby forming a masked data set for each source position, PA1 2) ramp filtering the cone beam projection data inside each masked data set, forming ramp filtered data sets, and PA1 3) subjecting the ramp filtered data sets to a weighted 3D backprojection into a 3D space corresponding to a complete field of view of a region of interest (ROI) of the object, thereby reconstructing in the 3D space a 3D image of the ROI in the object. PA1 m1) compute the integrals on line segments bounded by the mask, and PA1 m2) compute the difference between the line integrals so computed on adjacent parallel line segments.
The above prior art procedure will be referred to hereinafter as the 6-step process. It is assumed in this process that the entire cone beam image of the object is captured on the detector of the imaging system. Consider a plane Q(r,.theta.), which intersects the object, formed by the source and the line L(r,.theta.) on the detector at angle .theta. and at a distance r from the origin. Ignoring the function M(r,.theta.), the operations 1 through 6 compute the contribution to the reconstructed object density on the plane Q(r,.theta.) from the x-ray data illuminating the plane and its immediate vicinity. Since the 6-step process is detector driven, a contribution from the data illuminating the plane is computed every time the plane intersects the scan path and thus is illuminated by the x-ray beam. Consequently, the function M(r,.theta.) is used after the filter function in step 2 to normalize the results. Normalization is particularly undesirable since it requires pre-computing and storing a 2D array M(r,.theta.) for each source position along an imaging scan path. Since there are usually hundreds, if not thousands of source positions, this type of normalization is both computationally intensive and resource (computer memory) expensive. In the Kudo et al article, however, at page 203 it is noted that in the special case where the scan path is a circle, steps 1-5 can be simplified into a single convolution step, which essentially comprises ramp filtering the cone beam image in the direction of the scan path. This ramp filtering is equivalent to the well-known Feldkamp algorithm for a single circular orbit, such algorithm being described in particular detail in the article by L. A. Feldkamp, L. C. Davis, and J. W. Kress, entitled "Practical cone-beam algorithm" published in the J. Opt. Soc. Am. A. Vol. 1, 1984, pages 612-619, incorporated by reference herein (see in particular the convolution function equations 15 and 16 at page 614, which describe the convolution function as: ##EQU1## The key to this simplification is that in the special case of a circular scan path, the normalization function M(r,.theta.) is a constant, equal to 2. Consequently, the filtered projection at each r,.theta. that results after step 2, can merely be divided by 2 to compensate for the data redundancy.
In my prior U.S. patent application Ser. No 09/106,537 entitled SIMPLIFIED CONE BEAM IMAGE RECONSTRUCTION USING 3D BACKPROJECTION, filed Jun. 29, 1998 and incorporated herein by reference, I described how to use this ramp filter simplification for image reconstruction in a cone beam imaging system having source scan paths other than a single circle, such as a spiral scan path, and furthermore, in a cone beam imaging system having a short detector, i.e., a detector that captures only a portion of the cone beam image at each cone beam view (i.e., at each source position). As described in greater detail in my aforenoted prior U.S. Patent Application, my "3-step" simplified technique comprises:
A key part of my prior invention is the masking process. My prior invention allows the image reconstruction processing to be greatly speeded-up due to a reduced need for extensive computations, elimination of the requirement for normalization step 3 of the 6-step process (thereby obviating the need for significant memory allocation for normalization factors), and furthermore, the imaging system can use a short detector, i.e., one that does not acquire at each source position a complete view of the ROI of the object.
My prior invention can be considered as applying masking to steps 1 through 5 of the Kudo et al 6-step process, followed by step 6. The application of masking to Kudo et al's steps 1 and 2 is conceptually equivalent to the following operations:
It is intended that steps m1 and m2 yield a quantity proportional to the Radon derivative for the portion of plane Q(r,.theta.) defined by the current source position and the prior and the subsequent source positions. The angular ranges for various portions of plane Q(r,.theta.) are illustrated in FIG. 3, and will be described in greater detail later herein.
The operations in steps m1 and m2 are illustrated in FIG. 4. As shown therein, L, L.sub.1 ' and L.sub.2 ' are three closely spaced parallel line segments that are bound by a mask 400, where L is midway between L.sub.1 ' and L.sub.2 '. Line segment L is representative of many such line segments formed at various angles in mask 400, and corresponds to the previously described lines L (r,.theta.) of FIG. 1A and lines 8 of FIG. 1B, as well known to those skilled in this technology, which are used for computing Radon derivative data from the cone beam projection data. In the present technique, the arithmetic difference between the integrals computed for a given pair of line segments L.sub.1 ' and L.sub.2 ' in mask 400 is determined, and corresponds to the Radon derivative of the partial plane defined by the line segment L and the current source position, up to a multiplicative constant. In reality, however, this masking method yields only an approximation of the Radon derivative of the partial plane. This is because, as described in my recently issued U.S. Pat. No. 5,748,697, the Radon derivative for the relevant portion of plane Q(r,.theta.) should be computed as shown in FIG. 5 herein. As evident from FIG. 5, line segment L is the same as that shown in FIG. 4, however line segments L.sub.1 and L.sub.2 are obtained by orthogonal translation of L, as they should be. Accordingly, contrary to what is shown in FIG. 4, the ends of line segments L.sub.1 and L.sub.2 shown in FIG. 5 are not bound by mask 500. The ending of the line segments at the boundaries of the mask is referred to herein as "hard masking". As a consequence of the hard masking, the arithmetic difference between the integrals computed on line segments L'.sub.1 and L'.sub.2 does not yield an exact Radon derivative, i.e., one only differing only by a multiplicative constant.
It would be desirable to improve my simplified filtered backprojection image reconstruction processing technique so as to develop a more exact image reconstruction.